C   December, 2005 -- updated to work with MODFLOW-2005
C   September, 2000 -- updated to work with MODFLOW-2000
C   May, 2000 -- fixed error that caused incorrect critical head values
C   to be written in an external file when the option to write an
C   external file (IHCSV>0) is used.
C   June, 1996 -- 3 statements in the version documented
C   in TWRI 6-A2 have been modified in order to correct a problem.
C   Although subsidence is only meant to be active for layers in which
C   IBQ>0, some of the subroutines performed subsidence calculations when
C   IBQ<0.  Note that this was a problem only if negative IBQ values
C   were specified.  That is, the code has always worked correctly for
C   IBQ=0 and IBQ>0.
C   September, 2003 -- added the following:
C    1. Print a warning message that the IBS1 Package has been supseseded
C       by the SUB Package.
C    2. If the SUB Package and the IBS package are used simultaneously,
C       stop the simulation.
C
      MODULE GWFIBSMODULE
        INTEGER, SAVE, POINTER    ::IIBSCB,IIBSOC,ISUBFM,ICOMFM,IHCFM
        INTEGER, SAVE, POINTER    ::ISUBUN,ICOMUN,IHCUN
        INTEGER, SAVE,    DIMENSION(:),     POINTER ::IBQ
        INTEGER, SAVE,    DIMENSION(:),     POINTER ::IBQ1
        REAL,    SAVE,    DIMENSION(:,:,:), POINTER ::HC
        REAL,    SAVE,    DIMENSION(:,:,:), POINTER ::SCE
        REAL,    SAVE,    DIMENSION(:,:,:), POINTER ::SCV
        REAL,    SAVE,    DIMENSION(:,:,:), POINTER ::SUB
      TYPE GWFIBSTYPE
        INTEGER,POINTER    ::IIBSCB,IIBSOC,ISUBFM,ICOMFM,IHCFM
        INTEGER,POINTER    ::ISUBUN,ICOMUN,IHCUN
        INTEGER,    DIMENSION(:),     POINTER ::IBQ
        INTEGER,    DIMENSION(:),     POINTER ::IBQ1
        REAL,       DIMENSION(:,:,:), POINTER ::HC
        REAL,       DIMENSION(:,:,:), POINTER ::SCE
        REAL,       DIMENSION(:,:,:), POINTER ::SCV
        REAL,       DIMENSION(:,:,:), POINTER ::SUB
      END TYPE GWFIBSTYPE
      TYPE(GWFIBSTYPE),SAVE ::GWFIBSDAT(10)
      END MODULE GWFIBSMODULE



      SUBROUTINE GWF2IBS7AR(IN,INSUB,IGRID)
C
C  Based on
C-----VERSION 07JUN1996 GWF1IBS6ALP AND VERSION 1117 02JUN1988 GWF1IBS6RPP
C-----VERSION 01AUG1996 -- modified to allow 200 layers instead of 80
C     ******************************************************************
C     ALLOCATE ARRAY STORAGE FOR INTERBED STORAGE PACKAGE
C     READ INTERBED STORAGE DATA
C     ******************************************************************
C
C        SPECIFICATIONS:
C     ------------------------------------------------------------------
      USE GLOBAL,       ONLY: NCOL,NROW,NLAY,IOUT,DELR,DELC,HNEW
      USE GWFIBSMODULE, ONLY: HC,SCE,SCV,SUB,IBQ,IBQ1,IIBSCB,IIBSOC,
     1                        ISUBFM,ICOMFM,IHCFM,ISUBUN,ICOMUN,IHCUN
      CHARACTER*24 ANAME(4)
C
      DATA ANAME(1) /'   PRECONSOLIDATION HEAD'/
      DATA ANAME(2) /'ELASTIC INTERBED STORAGE'/
      DATA ANAME(3) /' VIRGIN INTERBED STORAGE'/
      DATA ANAME(4) /'     STARTING COMPACTION'/
C     ------------------------------------------------------------------
C
      ALLOCATE(IIBSCB,IIBSOC,ISUBFM,ICOMFM,IHCFM,ISUBUN,ICOMUN,IHCUN)
      ALLOCATE(IBQ(NLAY),IBQ1(NLAY))
C
C1------IDENTIFY PACKAGE.
      WRITE(IOUT,1)IN
    1 FORMAT(1H0,'IBS -- INTERBED STORAGE PACKAGE, VERSION 7,',
     1     ' 12/27/2005',' INPUT READ FROM UNIT',I3)
C1a------PRINT WARNING MESSAGE THAT PACKAGE HAS BEEN SUPERSEDED.
      WRITE(IOUT,103)
  103 FORMAT(1H0,'***NOTICE*** AS OF SEPTEMBER 2003, THE INTERBED ',
     1     'STORAGE PACKAGE HAS ',/,
     2     'BEEN SUPERSEDED BY THE SUBSIDENCE AND AQUIFER-SYSTEM ',
     3     'COMPACTION PACKAGE.',/,' SUPPORT FOR IBS MAY BE ',
     4     'DISCONTINUED IN THE FUTURE.')
C1b------PRINT A MESSAGE AND STOP THE SIMULATION OF BOTH IBS AND SUB
C1b------ ARE USED.
      IF(INSUB.GT.0) THEN
       WRITE(IOUT,104)
  104  FORMAT(1H0,'***ERROR*** THE IBS AND SUB PACKAGE SHOULD ',
     1     'NOT BOTH BE USED IN ',/,
     2     'THE SAME SIMULATION. ********STOPPING******** ')
       CALL USTOP(' ')
      ENDIF
C
C4------READ FLAG FOR STORING CELL-BY-CELL STORAGE CHANGES AND
C4------FLAG FOR PRINTING AND STORING COMPACTION, SUBSIDENCE, AND
C4------CRITICAL HEAD ARRAYS.
      READ(IN,3) IIBSCB,IIBSOC
    3 FORMAT(2I10)
C
C5------IF CELL-BY-CELL TERMS TO BE SAVED THEN PRINT UNIT NUMBER.
      IF(IIBSCB.GT.0) WRITE(IOUT,105) IIBSCB
  105 FORMAT(1X,'CELL-BY-CELL FLOW TERMS WILL BE SAVED ON UNIT',I3)
C
C5A-----IF OUTPUT CONTROL FOR PRINTING ARRAYS IS SELECTED PRINT MESSAGE.
      IF(IIBSOC.GT.0) WRITE(IOUT,106)
  106 FORMAT(1X,'OUTPUT CONTROL RECORDS FOR IBS PACKAGE WILL BE ',
     1 'READ EACH TIME STEP.')
C
C6------READ INDICATOR AND FIND OUT HOW MANY LAYERS HAVE INTERBED STORAGE.
      READ(IN,110) (IBQ(K),K=1,NLAY)
  110 FORMAT(40I2)
      NAQL=0
      DO 120 K=1,NLAY
      IF(IBQ(K).LE.0) GO TO 120
      NAQL=NAQL+1
      IBQ1(NAQL)=K
  120 CONTINUE
C
C7------IDENTIFY WHICH LAYERS HAVE INTERBED STORAGE.
      WRITE(IOUT,130) (IBQ1(K),K=1,NAQL)
  130 FORMAT(1X,'INTERBED STORAGE IN LAYER(S) ',80I2)
C
C8------ALLOCATE SPACE FOR THE ARRAYS HC, SCE, SCV, AND SUB.
      IF(NAQL.GT.0) THEN
        ALLOCATE(HC(NCOL,NROW,NAQL))
        ALLOCATE(SCE(NCOL,NROW,NAQL))
        ALLOCATE(SCV(NCOL,NROW,NAQL))
        ALLOCATE(SUB(NCOL,NROW,NAQL))
      ELSE
        ALLOCATE(HC(1,1,1))
        ALLOCATE(SCE(1,1,1))
        ALLOCATE(SCV(1,1,1))
        ALLOCATE(SUB(1,1,1))
      END IF
C
C9------READ IN STORAGE AND CRITICAL HEAD ARRAYS
      KQ=0
      DO 160 K=1,NLAY
      IF(IBQ(K).LE.0) GO TO 160
      KQ=KQ+1
      CALL U2DREL(HC(:,:,KQ),ANAME(1),NROW,NCOL,K,IN,IOUT)
      CALL U2DREL(SCE(:,:,KQ),ANAME(2),NROW,NCOL,K,IN,IOUT)
      CALL U2DREL(SCV(:,:,KQ),ANAME(3),NROW,NCOL,K,IN,IOUT)
      CALL U2DREL(SUB(:,:,KQ),ANAME(4),NROW,NCOL,K,IN,IOUT)
  160 CONTINUE
C
C10-----LOOP THROUGH ALL CELLS WITH INTERBED STORAGE.
      KQ=0
      DO 180 K=1,NLAY
      IF(IBQ(K).LE.0) GO TO 180
      KQ=KQ+1
      DO 170 IR=1,NROW
      DO 170 IC=1,NCOL
C
C11-----MULTIPLY STORAGE BY AREA TO GET STORAGE CAPACITY.
      AREA=DELR(IC)*DELC(IR)
      SCE(IC,IR,KQ)=SCE(IC,IR,KQ)*AREA
      SCV(IC,IR,KQ)=SCV(IC,IR,KQ)*AREA
C
C12-----MAKE SURE THAT PRECONSOLIDATION HEAD VALUES
C12-----ARE CONSISTANT WITH STARTING HEADS.
      IF(HC(IC,IR,KQ).GT.HNEW(IC,IR,K)) HC(IC,IR,KQ)=HNEW(IC,IR,K)
  170 CONTINUE
  180 CONTINUE
C
C13-----INITIALIZE AND READ OUTPUT FLAGS.
      ICOMFM=0
      ISUBFM=0
      IHCFM=0
      ICOMUN=0
      ISUBUN=0
      IHCUN=0
      IF(IIBSOC.LE.0) GO TO 200
      READ(IN,190) ISUBFM,ICOMFM,IHCFM,ISUBUN,ICOMUN,IHCUN
  190 FORMAT(6I10)
      WRITE(IOUT,191) ISUBFM,ICOMFM,IHCFM
  191 FORMAT(1H0,'    SUBSIDENCE PRINT FORMAT IS NUMBER',I4/
     1          '     COMPACTION PRINT FORMAT IS NUMBER',I4/
     2          '  CRITICAL HEAD PRINT FORMAT IS NUMBER',I4)
      IF(ISUBUN.GT.0) WRITE(IOUT,192) ISUBUN
  192 FORMAT(1H0,'    UNIT FOR SAVING SUBSIDENCE IS',I4)
      IF(ICOMUN.GT.0) WRITE(IOUT,193) ICOMUN
  193 FORMAT(1H ,'    UNIT FOR SAVING COMPACTION IS',I4)
      IF(IHCUN.GT.0)  WRITE(IOUT,194) IHCUN
  194 FORMAT(1H ,' UNIT FOR SAVING CRITICAL HEAD IS',I4)
C
C14-----RETURN
  200 CALL SGWF2IBS7PSV(IGRID)
      RETURN
      END
      SUBROUTINE GWF2IBS7ST(KPER,IGRID)
C
C-----Based on VERSION 15SEPT2000 GWF1IBS6ST
C     ******************************************************************
C     CHECK THAT NO STREE PERIOD IS STEADY STATE EXCEPT THE FIRST, AND
C     SET HC EQUAL TO THE STEADY-STATE HEAD IF STEADY-STATE HEAD IS
C     LOWER THAN HC.
C     ******************************************************************
C
C        SPECIFICATIONS:
C     ------------------------------------------------------------------
      USE GLOBAL,       ONLY: NCOL,NROW,NLAY,ISSFLG,HNEW,IOUT
      USE GWFIBSMODULE, ONLY:HC,IBQ
C     ------------------------------------------------------------------
C
C-------SET POINTERS FOR THE CURRENT GRID.
      CALL SGWF2IBS7PNT(IGRID)
C
C1------STOP if steady state after 1st stress period.
      IF(KPER.GT.1 .AND. ISSFLG(KPER).NE.0) THEN
        WRITE(IOUT,8)
    8   FORMAT(1X,'INTERBED STORAGE INAPPROPRIATE FOR A STEADY-STATE',
     1   /,1X,'STRESS PERIOD AFTER PERIOD 1 -- SIMULATION STOPPING.')
        CALL USTOP(' ')
      END IF
C
C4------MAKE SURE THAT PRECONSOLIDATION HEAD VALUES
C4------ARE CONSISTANT WITH STEADY-STATE HEADS.
      IF(KPER.EQ.2 .AND. ISSFLG(1).NE.0) THEN
        KQ=0
        DO 80 K=1,NLAY
        IF(IBQ(K).GT.0) THEN
          KQ=KQ+1
          DO 70 IR=1,NROW
          DO 70 IC=1,NCOL
          HTMP=HNEW(IC,IR,K)
          IF(HC(IC,IR,KQ).GT.HTMP) HC(IC,IR,KQ)=HTMP
   70     CONTINUE
        END IF
   80   CONTINUE
      END IF
C
      RETURN
      END
      SUBROUTINE GWF2IBS7FM(KPER,IGRID)
C
C  Based on:
C-----VERSION 07JUN1996 GWF1IBS6FM
C-----VERSION 01AUG1996 -- modified to allow 200 layers instead of 80
C     ******************************************************************
C        ADD INTERBED STORAGE TO RHS AND HCOF
C     ******************************************************************
C
C        SPECIFICATIONS:
C     ------------------------------------------------------------------
      USE GLOBAL,       ONLY:NCOL,NROW,NLAY,RHS,HCOF,HNEW,HOLD,IBOUND,
     1                       ISSFLG
      USE GWFBASMODULE, ONLY:DELT
      USE GWFIBSMODULE, ONLY:HC,SCE,SCV,IBQ
C     ------------------------------------------------------------------
C
C-------SET POINTERS FOR THE CURRENT GRID.
      CALL SGWF2IBS7PNT(IGRID)
C
C0------Return if stress period is steady state.
      ISSF=ISSFLG(KPER)
      IF(ISSF.NE.0) RETURN
C
C1------INITIALIZE
      TLED=1./DELT
      KQ=0
C
C2------FIND LAYERS WITH INTERBED STORAGE
      DO 110 K=1,NLAY
      IF(IBQ(K).LE.0) GO TO 110
      KQ=KQ+1
      DO 100 I=1,NROW
      DO 100 J=1,NCOL
      IF(IBOUND(J,I,K).LE.0) GO TO 100
C
C3------DETERMINE STORAGE CAPACITIES FOR CELL AT START AND END OF STEP
      RHO1=SCE(J,I,KQ)*TLED
      RHO2=RHO1
      HCTMP=HC(J,I,KQ)
      IF(HNEW(J,I,K).LT.HCTMP) RHO2=SCV(J,I,KQ)*TLED
C
C4------ADD APPROPRIATE TERMS TO RHS AND HCOF
      RHS(J,I,K)=RHS(J,I,K)-HCTMP*(RHO2-RHO1)-RHO1*HOLD(J,I,K)
      HCOF(J,I,K)=HCOF(J,I,K)-RHO2
  100 CONTINUE
  110 CONTINUE
C
C5------RETURN
      RETURN
      END
      SUBROUTINE GWF2IBS7BD(KSTP,KPER,IGRID)
C  Based on:
C-----VERSION 07JUN1996 GWF1IBS6BD
C-----VERSION 01AUG1996 -- modified to allow 200 layers instead of 80
C     ******************************************************************
C     CALCULATE VOLUMETRIC BUDGET FOR INTERBED STORAGE
C     ******************************************************************
C
C     SPECIFICATIONS:
C     ------------------------------------------------------------------
      USE GLOBAL,       ONLY:NCOL,NROW,NLAY,IBOUND,HOLD,HNEW,BUFF,
     1                       DELR,DELC,ISSFLG,IOUT
      USE GWFBASMODULE, ONLY:DELT,VBVL,VBNM,MSUM,ICBCFL
      USE GWFIBSMODULE, ONLY:HC,SCE,SCV,SUB,IBQ,IIBSCB
C
      CHARACTER*16 TEXT
      DATA TEXT /'INTERBED STORAGE'/
C     ------------------------------------------------------------------
C
C-------SET POINTERS FOR THE CURRENT GRID.
      CALL SGWF2IBS7PNT(IGRID)
      ISSF=ISSFLG(KPER)
C
C1------INITIALIZE CELL-BY-CELL FLOW TERM FLAG (IBD) AND
C1------ACCUMULATORS (STOIN AND STOUT).
      IBD=0
      STOIN=0.
      STOUT=0.
C
C2------TEST TO SEE IF CELL-BY-CELL FLOW TERMS ARE NEEDED.
      IF(ICBCFL.NE.0  .AND. IIBSCB.GT.0 ) THEN
C
C3------CELL-BY-CELL FLOW TERMS ARE NEEDED SET IBD AND CLEAR BUFFER.
        IBD=1
        DO 5 IL=1,NLAY
        DO 5 IR=1,NROW
        DO 5 IC=1,NCOL
        BUFF(IC,IR,IL)=0.
    5   CONTINUE
      END IF
C
C4------RUN THROUGH EVERY CELL IN THE GRID WITH INTERBED STORAGE.
      KQ=0
      TLED=1.
      IF(ISSF.EQ.0) TLED=1./DELT
      DO 110 K=1,NLAY
      IF(IBQ(K).LE.0 .OR. ISSF.NE.0) GO TO 110
      KQ=KQ+1
      DO 100 I=1,NROW
      DO 100 J=1,NCOL
C
C5------CALCULATE FLOW FROM STORAGE (VARIABLE HEAD CELLS ONLY)
      IF(IBOUND(J,I,K).LE.0) GO TO 100
      HHOLD=HOLD(J,I,K)
      HHNEW=HNEW(J,I,K)
      HHC=HC(J,I,KQ)
C
C6------GET STORAGE CAPACITIES AT BEGINNING AND END OF TIME STEP.
      SBGN=SCE(J,I,KQ)
      SEND=SBGN
      IF(HHNEW.LT.HHC) SEND=SCV(J,I,KQ)
C
C7------CALCULATE VOLUME CHANGE IN INTERBED STORAGE FOR TIME STEP.
      STRG=HHC*(SEND-SBGN)+SBGN*HHOLD-SEND*HHNEW
C
C8------ACCUMULATE SUBSIDENCE ASSOCIATED WITH CHANGE IN STORAGE
      SUB(J,I,KQ)=SUB(J,I,KQ)+STRG/(DELR(J)*DELC(I))
C
C9------IF C-B-C FLOW TERMS ARE TO BE SAVED THEN ADD RATE TO BUFFER.
      IF(IBD.EQ.1) BUFF(J,I,K)=BUFF(J,I,K)+STRG*TLED
C
C10-----SEE IF FLOW IS INTO OR OUT OF STORAGE.
      IF(STRG)94,100,96
   94 STOUT=STOUT-STRG
      GO TO 100
   96 STOIN=STOIN+STRG
  100 CONTINUE
  110 CONTINUE
C
C11-----IF C-B-C FLOW TERMS WILL BE SAVED CALL UBUDSV TO RECORD THEM.
      IF(IBD.EQ.1) CALL UBUDSV(KSTP,KPER,TEXT,IIBSCB,BUFF,NCOL,NROW,
     1                          NLAY,IOUT)
C
C12-----MOVE RATES,VOLUMES & LABELS INTO ARRAYS FOR PRINTING.
  200 VBVL(3,MSUM)=STOIN*TLED
      VBVL(4,MSUM)=STOUT*TLED
      VBVL(1,MSUM)=VBVL(1,MSUM)+STOIN
      VBVL(2,MSUM)=VBVL(2,MSUM)+STOUT
      VBNM(MSUM)=TEXT
C
C13-----INCREMENT BUDGET TERM COUNTER
      MSUM=MSUM+1
      IF(ISSF.NE.0) RETURN
C
C14-----UPDATE PRECONSOLIDATION HEAD ARRAY
      KQ=0
      DO 310 K=1,NLAY
      IF(IBQ(K).LE.0) GO TO 310
      KQ=KQ+1
      DO 300 I=1,NROW
      DO 300 J=1,NCOL
      IF(IBOUND(J,I,K).LE.0) GO TO 300
      HHNEW=HNEW(J,I,K)
      IF(HHNEW.LT.HC(J,I,KQ)) HC(J,I,KQ)=HHNEW
  300 CONTINUE
  310 CONTINUE
C
C15-----RETURN
      RETURN
      END
      SUBROUTINE GWF2IBS7OT(KSTP,KPER,IN,IGRID)
C  Based on
C-----VERSION 07JUN1996 GWF1IBS6OT
C-----VERSION 01AUG1996 -- modified to allow 200 layers instead of 80
C     ******************************************************************
C     PRINT AND STORE SUBSIDENCE, COMPACTION AND CRITICAL HEAD.
C     ******************************************************************
C
C     SPECIFICATIONS:
C     ------------------------------------------------------------------
      USE GLOBAL,       ONLY:NCOL,NROW,NLAY,BUFF,NSTP,ISSFLG,IOUT
      USE GWFBASMODULE, ONLY:PERTIM,TOTIM
      USE GWFIBSMODULE, ONLY:IBQ,SUB,HC,IIBSOC,ISUBFM,ICOMFM,IHCFM,
     1                          ISUBUN,ICOMUN,IHCUN
      CHARACTER*16 TEXT(3)
      DATA TEXT(1) /'      SUBSIDENCE'/
      DATA TEXT(2) /'      COMPACTION'/
      DATA TEXT(3) /'   CRITICAL HEAD'/
C     ------------------------------------------------------------------
C
C-------SET POINTERS FOR THE CURRENT GRID.
      CALL SGWF2IBS7PNT(IGRID)
      ISSF=ISSFLG(KPER)
      IF(ISSF.NE.0) RETURN
C
C1------INITIALIZE FLAGS FOR PRINTING AND SAVING SUBSIDENCE, COMPACTION,
C1------AND CRITICAL HEAD
      ISUBPR=0
      ICOMPR=0
      IHCPR=0
      ISUBSV=0
      ICOMSV=0
      IHCSV=0
      IF(KSTP.EQ.NSTP(KPER)) ISUBPR=1
C2------READ FLAGS FOR PRINTING AND SAVING.
      IF(IIBSOC.LE.0) GO TO 28
      READ(IN,10) ISUBPR,ICOMPR,IHCPR,ISUBSV,ICOMSV,IHCSV
   10 FORMAT(6I10)
      WRITE(IOUT,15) ISUBPR,ICOMPR,IHCPR,ISUBSV,ICOMSV,IHCSV
   15 FORMAT(1H0,'FLAGS FOR PRINTING AND STORING SUBSIDENCE, ',
     1 'COMPACTION, AND CRITICAL HEAD:'/
     2 '   ISUBPR    ICOMPR    IHCPR     ISUBSV    ICOMSV    IHCSV   '/
     3 ' ------------------------------------------------------------'/
     4 I6,5I10)
C
C3------PRINT AND STORE SUBSIDENCE, FIRST, CLEAR OUT BUFF.
   28 IF(ISUBPR.LE.0.AND.ISUBSV.LE.0) GO TO 100
      DO 30 IR=1,NROW
      DO 30 IC=1,NCOL
      BUFF(IC,IR,1)=0.
   30 CONTINUE
C
C4------SUM COMPACTION IN ALL LAYERS TO GET SUBSIDENCE.
      KQ=0
      DO 50 K=1,NLAY
      IF(IBQ(K).LE.0) GO TO 50
      KQ=KQ+1
      DO 40 I=1,NROW
      DO 40 J=1,NCOL
      BUFF(J,I,1)=BUFF(J,I,1)+SUB(J,I,KQ)
   40 CONTINUE
   50 CONTINUE
C
C5------PRINT SUBSIDENCE.
      IF(ISUBPR.LE.0) GO TO 60
      IF(ISUBFM.LT.0) CALL ULAPRS(BUFF,TEXT(1),KSTP,KPER,NCOL,NROW,1,
     1          -ISUBFM,IOUT)
      IF(ISUBFM.GE.0) CALL ULAPRW(BUFF,TEXT(1),KSTP,KPER,NCOL,NROW,1,
     1           ISUBFM,IOUT)
C
C6------STORE SUBSIDENCE.
   60 IF(ISUBSV.LE.0) GO TO 100
      CALL ULASAV(BUFF,TEXT(1),KSTP,KPER,PERTIM,TOTIM,NCOL,NROW,1,
     1             ISUBUN)
C
C7------PRINT COMPACTION FOR ALL LAYERS WITH INTERBED STORAGE.
  100 IF(ICOMPR.LE.0) GO TO 140
      KQ=0
      DO 130 K=1,NLAY
      IF(IBQ(K).LE.0) GO TO 130
      KQ=KQ+1
      IF(ICOMFM.LT.0) CALL ULAPRS(SUB(:,:,KQ),TEXT(2),KSTP,KPER,NCOL,
     1          NROW,K,-ICOMFM,IOUT)
      IF(ICOMFM.GE.0) CALL ULAPRW(SUB(:,:,KQ),TEXT(2),KSTP,KPER,NCOL,
     1           NROW,K,ICOMFM,IOUT)
  130 CONTINUE
C
C8------SAVE COMPACTION FOR ALL LAYERS WITH INTERBED STORAGE.
  140 IF(ICOMSV.LE.0) GO TO 200
      KQ=0
      DO 160 K=1,NLAY
      IF(IBQ(K).LE.0) GO TO 160
      KQ=KQ+1
      CALL ULASAV(SUB(:,:,KQ),TEXT(2),KSTP,KPER,PERTIM,TOTIM,NCOL,
     1            NROW,K,ICOMUN)
  160 CONTINUE
C
C9------PRINT CRITICAL HEAD FOR ALL LAYERS WITH INTERBED STORAGE.
  200 IF(IHCPR.LE.0) GO TO 240
      KQ=0
      DO 230 K=1,NLAY
      IF(IBQ(K).LE.0) GO TO 230
      KQ=KQ+1
      IF(IHCFM.LT.0) CALL ULAPRS(HC(:,:,KQ),TEXT(3),KSTP,KPER,NCOL,
     1          NROW,K,-IHCFM,IOUT)
      IF(IHCFM.GE.0) CALL ULAPRW(HC(:,:,KQ),TEXT(3),KSTP,KPER,NCOL,
     1           NROW,K,IHCFM,IOUT)
  230 CONTINUE
C
C10-----SAVE CRITICAL HEAD FOR ALL LAYERS WITH INTERBED STORAGE.
  240 IF(IHCSV.LE.0) GO TO 300
      KQ=0
      DO 260 K=1,NLAY
      IF(IBQ(K).LE.0) GO TO 260
      KQ=KQ+1
      CALL ULASAV(HC(:,:,KQ),TEXT(3),KSTP,KPER,PERTIM,TOTIM,NCOL,
     1            NROW,K,IHCUN)
  260 CONTINUE
C
C11-----RETURN
  300 RETURN
      END
      SUBROUTINE GWF2IBS7DA(IGRID)
C  Deallocate IBS DATA
      USE GWFIBSMODULE
C
        DEALLOCATE(GWFIBSDAT(IGRID)%IIBSCB)
        DEALLOCATE(GWFIBSDAT(IGRID)%IIBSOC)
        DEALLOCATE(GWFIBSDAT(IGRID)%ISUBFM)
        DEALLOCATE(GWFIBSDAT(IGRID)%ICOMFM)
        DEALLOCATE(GWFIBSDAT(IGRID)%IHCFM)
        DEALLOCATE(GWFIBSDAT(IGRID)%ISUBUN)
        DEALLOCATE(GWFIBSDAT(IGRID)%ICOMUN)
        DEALLOCATE(GWFIBSDAT(IGRID)%IHCUN)
        DEALLOCATE(GWFIBSDAT(IGRID)%IBQ)
        DEALLOCATE(GWFIBSDAT(IGRID)%IBQ1)
        DEALLOCATE(GWFIBSDAT(IGRID)%HC)
        DEALLOCATE(GWFIBSDAT(IGRID)%SCE)
        DEALLOCATE(GWFIBSDAT(IGRID)%SCV)
        DEALLOCATE(GWFIBSDAT(IGRID)%SUB)
C
      RETURN
      END
      SUBROUTINE SGWF2IBS7PNT(IGRID)
C  Set IBS pointers for grid.
      USE GWFIBSMODULE
C
        IIBSCB=>GWFIBSDAT(IGRID)%IIBSCB
        IIBSOC=>GWFIBSDAT(IGRID)%IIBSOC
        ISUBFM=>GWFIBSDAT(IGRID)%ISUBFM
        ICOMFM=>GWFIBSDAT(IGRID)%ICOMFM
        IHCFM=>GWFIBSDAT(IGRID)%IHCFM
        ISUBUN=>GWFIBSDAT(IGRID)%ISUBUN
        ICOMUN=>GWFIBSDAT(IGRID)%ICOMUN
        IHCUN=>GWFIBSDAT(IGRID)%IHCUN
        IBQ=>GWFIBSDAT(IGRID)%IBQ
        IBQ1=>GWFIBSDAT(IGRID)%IBQ1
        HC=>GWFIBSDAT(IGRID)%HC
        SCE=>GWFIBSDAT(IGRID)%SCE
        SCV=>GWFIBSDAT(IGRID)%SCV
        SUB=>GWFIBSDAT(IGRID)%SUB
C
      RETURN
      END
      SUBROUTINE SGWF2IBS7PSV(IGRID)
C  Save IBS pointers for grid.
      USE GWFIBSMODULE
C
        GWFIBSDAT(IGRID)%IIBSCB=>IIBSCB
        GWFIBSDAT(IGRID)%IIBSOC=>IIBSOC
        GWFIBSDAT(IGRID)%ISUBFM=>ISUBFM
        GWFIBSDAT(IGRID)%ICOMFM=>ICOMFM
        GWFIBSDAT(IGRID)%IHCFM=>IHCFM
        GWFIBSDAT(IGRID)%ISUBUN=>ISUBUN
        GWFIBSDAT(IGRID)%ICOMUN=>ICOMUN
        GWFIBSDAT(IGRID)%IHCUN=>IHCUN
        GWFIBSDAT(IGRID)%IBQ=>IBQ
        GWFIBSDAT(IGRID)%IBQ1=>IBQ1
        GWFIBSDAT(IGRID)%HC=>HC
        GWFIBSDAT(IGRID)%SCE=>SCE
        GWFIBSDAT(IGRID)%SCV=>SCV
        GWFIBSDAT(IGRID)%SUB=>SUB
C
      RETURN
      END
